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5.2 Accuracy and Stability for Ut = cu x 

This section begins a major topic in scientific computing: Initial-value problems 
for partial differential equations. Naturally we start with linear equations that 
involve only one space dimension x (and time t). The exact solution is u(x, t) and its 
discrete approximation on a space-time grid has the form Uj )n = U(jAx,nAt). We 
want to know if U is near u — how close they are and how stable U is. 

Begin with the simplest wave equation (first-order, linear, constant coefficient): 

du Ou 

One-way wave equation — = c— . (1) 

at ox 

We are given u(x, 0) at time t — 0. We want to find u(x, t) for all t > 0. For simplicity, 
these functions are defined on the whole line — oo < x < oo. There are no difficulties 
with boundaries (where waves could change direction and bounce back). 

The solution u(x,t) will have the typical feature of hyperbolic equations: signals 
travel at finite speed. Unlike the second-order wave equation u u = c 2 u xx , this first- 
order equation u t = cu x sends signals in one direction only. 



Solution for u(x,0) = e lkx 

Throughout this chapter I will solve for a pure exponential u(x, 0) = e tkx . At every 
time t, the solution remains a multiple Ge tkx . The growth factor G will depend 
on the frequency k and the time t, but different frequencies do not mix. Substituting 
u = G(k,t)e lkx into u t = cu x yields a simple ordinary differential equation for G, 
because we can cancel e lkx . The derivative of e lkx produces the factor ik: 

dG dG 
u t = cu x is —e ikx = ikcGe ikx or — = ikcG . (2) 
dt dt 

The growth factor is G(k, t) = e zkct . The initial value is G = 1. 

An exponential solution to — = c— is u(x, t) = e ikct e ikx = e ik i x + c£ ) . (3) 

ut ox 

Immediately we see two important features of this solution: 

1. The growth factor G = e lkct has absolute value |G| = 1. 

2. The initial function e tkx moves to the left with fixed velocity c, to e lk ^ x+ct \ 

The initial value at the origin is w(0,0) = e lk0 = 1. This value u = 1 appears 
at all points on the line x + ct = 0. The initial data propagates along the 
characteristic lines x + ct = constant, in Figure 5.3. Right now we know this 
fact for the special solutions e lk ( x+ct )_ Soon we will know it for all solutions. 
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Figure 5.3: The solution u(x,t) moves left with speed c, along characteristic lines. 



Figure 5.3 shows the travel path of the solution in the x-t plane (we are introducing 
the characteristic lines). Figure 5.4 will graph the solution itself at times and t. 
That step function combines exponentials e lk ( x+ct ) for different frequencies k. By 
linearity we can add those solutions. 



Solution for Every u(x,0) 

In almost all partial differential equations, the solution changes shape as it travels. 
Here the shape stays the same. All pure exponentials travel at the same velocity c, 
so every initial function moves with that velocity. We can write down the solution: 

On du 

General solution — = c — is solved by u(x, t) = u(x + ct, 0) . (4) 

The solution is a function only of x + ct. That makes it constant along characteristic 
lines, where x + ct is constant. This dependence on x + ct also makes it satisfy the 
equation m = cu x , by the chain rule. If we take u = (x + ct) n as an example, the 
extra factor c appears in du/dt: 

du . , , du . .„_, , . , . du 

— = n(x + ct) and — = cn(x + ct) which is c— . 
ox at ox 

A Taylor series person would combine those powers (different n) to produce a large 
family of solutions. A Fourier series person combines exponentials (different k) to 
produce an even larger family. In fact all solutions are functions of x + ct alone. 

Here are two important initial functions — a light flashes or a dam breaks. 

Example 1 u(x,0) = delta function 5(x) = flash of light at x = 0,t = 

By our formula (4), the solution is u(x, t) = 8(x + ct). The light flash reaches the point 
x = — c at the time t = 1. It reaches x = —2c at the time t = 2. The impulse is traveling 
to the left at speed \dx/dt\ = c. In this example all frequencies k are present in equal 
amounts, because the Fourier transform of a delta function is a constant. 
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Notice that a point goes dark again as soon as the flash passes through. This is the 
Huygens principle in 1 and 3 dimensions. If we lived in two or four dimensions, the wave 
would not pass all at once and we wouldn't see clearly. 



Example 2 u(x,0) = step function S(x) = wall of water at x = 0,t = 

The solution S(x + ct) is the moving step function in Figure 5.4. The wall of water travels 
to the left (one-way wave). At time t, the "tsunami" reaches the point x = —ct. The 
flash of light will get there first, because its speed c is greater than the tsunami speed. 
That is why a warning is possible for an approaching tsunami. 



u(x, 0) 



u(x, t) 







x 



x 



-ct 
later profile S(x + ct) at time t 



initial profile S(x) at t — 
Figure 5.4: The wall travels left with velocity c (all waves e lkx do too). 



An actual tsunami is described by the nonlinear "shallow water equations" that come 
later. The feature of finite speed still holds. 



Finite Difference Methods for Ut = cu x 

The one-way wave equation is a perfect example for creating and testing finite dif- 
ference approximations. We can replace du/dt by a forward difference with step At. 
Here are four choices for the discrete form of du/dx at meshpoint iAx: 

1. Forward = — — = upwind: Low accuracy, conditionally stable for c > 0. 

Ax 

2. Centered = — 1 : Unstable after a few steps as we will prove! 

3. Lax-Friedrichs: (20) has low accuracy, conditionally stable also for c < 0. 

4. Lax- Wendroff : (14) has extra accuracy, conditionally stable also for c < 0. 

The list doesn't end there. We have reached a central problem of scientific computing, 
to construct approximations that are stable and accurate and fast. That topic can't 
be developed on one page, especially when we move to nonlinear equations. 

Conditionally stable means that the time step At is restricted. The need for this 
restriction was noticed by Courant, Friedrichs, and Lewy. When the space difference 
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reaches no further than x + Ax, there is an automatic stability restriction: 

At 

CFL requirement for stability r = c < 1 . (5) 

Ax 

That number cAt/Ax is often called the Courant number. (It was really Lewy who 
recognized that r < 1 is necessary for stability and convergence.) The reasoning is 
straightforward, based on using the initial value that controls u(x,t): 

The true solution at (x,t) equals the initial value u(x + ct,0). Taking n 
discrete steps to reach t = n At uses information on the initial values as far 
out as x + n Ax. If x + ct is further than x + n Ax, the method can't work: 

CFL condition x + ct<x + nAx or cnAt<nAx or r = c— — <1. (6) 

If the difference equation uses U(x + 2 Ax, t), then CFL relaxes to r < 2. 

A particular finite difference equation might require a tighter restriction on At for 
stability It might even be unstable for all ratios r (we hope not). The only route to 
unconditional stability for all At is an implicit method, which computes x-differences 
at the new time t + At. This will be useful later for diffusion terms like u xx . For 
advection terms (first derivatives), explicit methods with a CFL limitation are usually 
accepted because a much larger At would lose accuracy as well as stability. 

To repeat, if r > 1 then the finite difference solution at x,t does not use initial 
value information near the correct point x* = x + ct. Hopeless. 



Accuracy of the Upwind Difference Equation 

Linear problems with constant coefficients are the ones to understand first. Exactly 
as for differential equations, we can follow each pure exponential e tkx . After a single 
time step, there will be a growth factor in U(x,At) = Ge lkx . That growth factor 
G(k, At, Ax) may have magnitude \G\ < 1 or \G\ > 1. This will control stability or 
instability. The order of accuracy (if we compute in the k-u domain) comes from 
comparing G with the true factor e lkcAt from the differential equation. 

We now determine that the order of accuracy is p = 1 for the upwind method. 

Forward differences m*,t + *) - U(*,t) =c U(* + **,t) -Ufrt) 

At Ax w 

We will test the accuracy in the x-t domain and then the k-u domain. Either way we 
use Taylor series to check the leading terms. Substituting the true solution u(x, t) in 
place of U(x,t), its forward differences are 

Time ^- [u(x, t + At) - u(x, t)) = u t + -Atut t -\ (8) 

c 1 
Space ^— [u(x + Ax, t) — u(x, t)] = cu x + -c Axu xx + • • • (9) 

On the right side, u t = cu x is good. One more derivative gives u tt = cu xt = c 2 u xx . 
Notice c 2 . Then Atuu matches cAxu xx only in the special case cAt = Ax: 
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1 9 1 cAt 

-Ak M ra equals -cAxw^ only if r = = 1. 

For any ratio r ^ 1, the difference between (8) and (9) has a first-order error. Let 
me show this also in the k-u Fourier picture and then improve to second-order. 

Fix the ratio r = cAt/Ax as Ax — > and At — > 0. In the difference equation (7), 
write each new value at time t + At as a combination of two old values of U: 

Difference equation U(x, t + At) = (1 — r) U(x, t)+r U(x + Ax, t) . (10) 

Starting from U(x, 0) = e tkx we quickly find the growth factor G at time At: 

After 1 step (1 - r)e ikx + r e ik ^ x+Ax) = [l - r + r e ifeAx ] e ikx = Ge ikx . (11) 

To test the accuracy, compare this G = G approx to the exact growth factor e ichAt . 
Use the power series 1 + x + x 2 /2! + • • • for any e x : 

AccuracyG approx = 1 - r + r e tkAx = (1 — r) + r + r(ikAx) + (ikAx) 2 H 

Gexact = e icfcA * = e irfcAx = 1 + irk Ax + i(zr£;A£) 2 + ■ • • (12) 

The first terms agree as expected. Forward differences replaced derivatives, and the 
method is consistent. We saw Ut = cu x in comparing (8) with (9). The next terms 
do not agree unless r = r 2 : 

Compare ^r(ikAx) 2 with ^r 2 (ikAx) 2 . Single-step error of order (kAt) 2 . (13) 

After 1/At steps, those errors of order k 2 (At) 2 give a final error 0{k 2 At). Forward 
differences are only first order accurate, and so is the whole method. 

The special case r = 1 means cAt = Ax. The difference equation is exactly 
correct. The true and approximate solutions at (x, At) are both u(x + Ax,0). We 
are on the characteristic line in Figure 5.3. This is an interesting special case (the 
golden At, but hard to repeat in scientific computing when c varies). 

Conclusion Except when r = r 2 , the upwind method is first-order accurate. 



Higher Accuracy for Lax-Wendroff 

To upgrade the accuracy, we need to match the | At u u error term in the forward time 
difference by an additional space difference that gives \Atc 2 u xx . This is achieved by 
the Lax-Wendroff method: 

U(x, t + At) - U(x, t) _ U(x + Ax, t) - U(x - Ax, t) 

At ~ ° 2Ar . .v 

At 2 (U{x + Ax,t) -2U(x,t) + U(x- Ax,t)^ 1 ' 

+ T° V (Ax) 2 

Substituting the true solution, that second difference produces \c 2 Atu xx plus higher 
order terms. This cancels the \Atuu error term in the time difference, computed 
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in equation (8). (Remember u u = cu xt = c 2 u xx . The centered difference has no Ax 
term.) Thus Lax-Wendroff has second-order accuracy. 

To see this in the k-uj frequency domain, rewrite the LW difference equation (14): 

U(x,t + At) = {l-r 2 )U{x,t) + ^{r 2 + r)U(x + Ax,t) + ^(r 2 -r)U{x-Ax,t). (15) 

Substitute U(x, t) = e lkx to find the one-step growth factor G at time t + At: 

Growth factor for LW G = (1 - r 2 ) + -(r 2 + r)e ikAx + -(r 2 - r)e~ ikAx . 

2 2 

Expanding e lkAx and e~ ikAx in powers of ik Ax, this becomes 

G = 1 + r(ikAx) + ^r 2 (ikAx) 2 + 0(kAxf . (16) 

Comparing with G e xact = e %rkAx in equation (12), three terms agree. So the one-step 
error is of order (kAx) 3 . After 1/At steps the second-order accuracy of Lax-Wendroff 
is confirmed. 

Figure 5.5 shows by actual computation the improvement in accuracy. For a first- 
order method, the "wall of water" is smeared out. High frequencies have growth 
factors | much smaller than 1. There is too much dissipation. For the first-order 
Lax-Friedrichs method, the dissipation is even worse (Problem 2). The second-order 
Lax-Wendroff method stays much closer to the discontinuity. But it's not perfect— 
those oscillations are not good. 

For an ideal difference equation, we want to add enough dissipation very close to 
the shock, to avoid that oscillation (the Gibbs phenomenon). A lot of thought has 
gone into high resolution methods, to capture shock waves cleanly. 

Greater accuracy is achievable by including more terms in the difference equation. 
If we go from the three terms in Lax-Wendroff to five terms, we can reach fourth- 
order accuracy. If we use all values U (j Ax, nAt) at every time step, which requires 
more work, we can achieve spectral accuracy. Then the error decreases faster than 
any power of Ax, provided u(x, t) is smooth enough to allow derivatives of all orders. 
Section gives a separate discussion of this spectral method. 



Stability of the Four Finite Difference Methods 

Now we turn from accuracy to stability. Accuracy requires G to stay close to the 
true e tckAt . Stability requires G to stay inside the unit circle. We need \G\ < 1 for all 
frequencies k or the finite difference approximation G n e lkx will blow up. 

We now check whether or not |G| < 1, in the four methods. 



1. Forward differences in space and time: At// At = cAU/Ax (upwind). 
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Recall from equation (11) that G = 1 — r + re lkAx . If the Courant number r is 
between and 1, the triangle inequality gives \G\ < 1: 

Stability for < r < 1 |G| < |1 - r| + \re ikAx \ = 1 - r + r = 1 . (17) 

This sufficient condition < cAt/Ax < 1 is exactly the same as the Courant- 
Friedrichs-Lewy necessary condition! They reasoned that U(x,nAt) depends on the 
initial values between x and x + nAx. That domain of dependence must include 
the point x + cnAt. (Otherwise, changing the initial value at the point x + cnAt 
would change the true solution u but not the approximation U .) Then cnAt must 
lie between and nAx, which means that < r < 1. 

Figure 5.5 shows G in the stable case r = | and the unstable case r = | (when At 
is too large). As k varies, and e tkAx goes around a unit circle, the complex number 
G = 1 — r + re tkAx goes in a circle of radius r. The center is 1 — r. Always G — 1 at 
zero frequency (constant solution, no growth). 




Figure 5.5: Stable (upwind) and unstable (centered): CFL numbers r = | and r = |. 



2. Forward difference in time, centered difference in space. 

This combination is never stable! The shorthand Uj n will stand for U(jAx,nAt): 

Uj,n+1 — Uj,n _ Uj+l,n ~ Uj-l,n TJ _ TJ , T (tt tt \ (-io\ 

At ~ 2Ax u i,n+i - Uj.n + ^ K u j+i,n ~ Uj-\,n) ■ {^) 

Those coefficients 1 and r/2 and — r/2 go into the growth factor G, when the solution 
is a pure exponential and e lkx is factored out: 

Unstable: \G\ > 1 G = 1 + ^e ikAx - L e ~ ikAx = 1 + ir sin fcA^ . (19) 

The real part is 1. The magnitude is \G\ > 1. Its graph is on the left side of Figure 5.6. 
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3. Lax-Friedrichs Method (upwind-downwind). 

We can recover stability for centered differences by changing the time difference. 
Replace Uj >n by the average |(Z7j + i )n + Z7j_i jn ) of its neighbors: 

Lax-Friedrichs U hn+l - \{U +l , n + U^ n ) = U s+1 , n - U^, n 

At 2Ax K ' 

Two old values Uj + i jU and £/j-i, n produce each new value Uj >n +i. Moving terms to 
the right-hand side, the coefficients are |(1 + r) and |(1 — r). The growth factor is 

G = e ikAx + ^—^ e~ ikAx = cos kAx + ir sin kAx . (21) 

The absolute value is \G\ 2 = (coskAx) 2 + r 2 (sin k Ax) 2 . In Figure 5.6, \G\ < 1 when 
r 2 < 1. This stability condition agrees again with the CFL condition. 



Forward in time— centered in space 

<— G = 1 + ir sin kAx 



\G\>1 
Unstable 



Lax-Friedrichs 

G = cos kAx + ir sin kAx 
Stable for r 2 < 1 



G|<1, 



Figure 5.6: Equation (18) is unstable for all r. Equation (20) is stable for r 2 < 1. 

Notice that c and r can be negative. The wave can go either way! This will 
be useful for the two-way wave equation, but the accuracy is still first-order. The 
Lax-Friedrichs G matches the next term in the exact growth factor only if r 2 = 1: 



G = cos kAx + ir sin kAx = 1 + irk Ax — -(kAx) 2 + 



(22) 



G 



exact = e ikrAx = l + irkAx + -i 2 r 2 {kAxf + 



In the exceptional cases r = 1 and r = — 1, G agrees with G ex act- Staying exactly 
on the characteristic line, t/^n+i matches the true u(jAx,t + At). For r 2 < 1, Lax- 
Friedrichs has an important advantage and disadvantage: 

Good Each new C/j, n +i is a positive combination of old values. 
Not good The accuracy is only first-order. 



Problem 6 will show that second-order is impossible with positive coefficients. 
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4. Lax-Wendroff Method (second-order accurate). 

The LW difference equation (14) combines Uj >n and E/j-i,n an d ^+i,n to compute the 
new value Uj >n+ i. The coefficients of these three old values go into G: 

Lax-Wendroff G = (1 — r 2 ) + -(r 2 + r)e ikAx + -(r 2 — r)e~ ikAx . (23) 

This is G = 1 — r 2 + r 2 cos A;Ax + sin fcAx. At the dangerous frequency kAx = tt, 
the growth factor is 1 — 2r 2 . That stays above —1 if r 2 < 1. 

Problem 5 shows that |G| < 1 for every kAx. Lax-Wendroff is stable when- 
ever the CFL condition r 2 < 1 is satisfied. Again the wave can go either way 
(or both ways) since c and r can be negative. This is the most accurate of the five 
methods in Figure 5.7. 



upwind 
stable 

if r < 1 



wrong way 
unstable 
all At 



centered 
unstable 
all At 



Lax-Friedrichs 
stable 

if \r\ < 1 



Lax-Wendroff 
stable 

if Irl < 1 



Figure 5.7: Difference methods for the one-way wave equation u% = cu a 



Equivalence of Stability and Convergence 

Does the discrete solution U approach the true solution u as At — > ? The ex- 
pected answer is yes. But there are two requirements for convergence, and one of 
them — stability — is by no means automatic. The other requirement is consistency— 
the discrete problem must approximate the correct continuous problem. The fact that 
these two properties are sufficient for convergence, and also necessary for convergence, 
is the fundamental theorem of numerical analysis: 

. . L . Stability is equivalent to convergence, for a consistent 

Lax equivalence theorem . . ,, ... 

approximation to a well-posed linear problem. 

Lax proved the equivalence theorem for initial- value problems. The rate of conver- 
gence is given in (26). The theorem is equally true for boundary- value problems, and 
for the approximation of functions, and for the approximation of integrals. It applies 
to every discretization, when the given problem Lu = f is replaced by L h U h = f h . 
Assuming the inputs / and fh are close, we will prove that u and Uh are close — 
provided Lh is stable. The key points of the proof take only a few lines when the 
equation is linear, and you will see the essence of this fundamental theorem. 

Suppose / is changed to fh and L is replaced by Lh- The requirements are 
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Consistency: fh~* f an d LhU — > Lu for smooth solutions u. 
Well-posed: The inverse of L is bounded: \\u\\ = — C||/ll- 

Stability: The inverses L^ 1 remain uniformly bounded: HL^ 1 /?*!! < C||/^||. 

Under those conditions, the approximation Uh = L^ fh will approach u as h goes to 
zero. We subtract and add L^Lu = L h l f when u is smooth: 

Convergence u-U h = L h x (L h u - Lu) + L h 1 (/ - f h ) -> . (24) 

Consistency controls the quantities in parentheses (they go to zero). Stability controls 
the operators L^ 1 that act on them. Well-posedness controls the approximation of 
all solutions by smooth solutions. Then always Uh converges to u. 

If stability fails, there will be an input for which the approximations Uh = Lj^f 
are not bounded. The uniform boundedness theorem produces this bad f, from the 
inputs fh on which instability gives H-Z^" //Jl — > oo. Convergence fails for this /. 

A perfect equivalence theorem goes a little further, after careful definitions: 

Consistency + Stability -<=>- Well-posedness + Convergence . 

Our effort will now concentrate on initial-value problems, to estimate the error (the 
convergence rate) in u — Uh- The parameter h becomes At. We take n steps. 

The Rate of Convergence 

Consistency means that the error at each time step goes to zero as the mesh is refined. 
Our Taylor series estimates have done more: The order of accuracy gives the rate 
that this one-step error goes to zero. The problem is to extend this local rate to a 
global rate of convergence, accumulating the errors over n time steps. 

Let me write S for a single finite difference step, so U(t + At) = S U(t). The 
corresponding step for the differential equation will be u(t + At) = Ru(t). Then 
consistency means that Su is close to Ru, and the order of accuracy p tells how close: 

Accuracy of discretization \\Su — Ru\\ < Ci(At) p+1 for smooth solutions u. 
Well-posed problem ^ C 2 ||tfc|| for n At < T. 

Stable approximations || S n U\\ < C 3 \\U\\ for nAt<T. 

The difference between U = S n u(0) and the true u = R n u( ) is (S n — R n )u( ). 
The key idea is a "telescoping identity" that involves n single-step differences S — R: 

S n - R n = S n -\S -R) + S n ~ 2 {S - R)R + ■ ■ ■ + {S - R)R n - 1 . (25) 

Each of those n terms has a clear meaning. First, a power R k carries u(0) to the true 
solution u(kAt). Then (S — R)u(k At) gives the error at step k of order (At) p+1 . 
Then powers of S carry that one-step error forward to time nAt. By stability, this 
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amplifies the error by no more than C3. There are n <T/ At steps. The final rate 
of convergence for smooth solutions is (At) p : 



\\U(nAt) - u{nAt)\\ = \\(S n - R n )u(0)\\ < CiC 2 C 3 ^(A*) P+1 = CiC 2 C 3 T(At)* . 

(26) 



Notice how smoothness was needed in the Taylor series (8) and (9), when At and 
Aa; multiplied u u and u xx . That first-order accuracy would not apply if u or u t or 
u x had a jump. Still the order of accuracy p gives a practical estimate of the overall 
approximation error u — U . The problem of scientific computing is to get beyond 
p = 1 while maintaining stability and speed. 



1 Integrate u t = cu x from —00 to 00 to prove that mass is conserved: dM/dt = 0. 
Multiply by u and integrate uu t = c uu x to prove that energy is also conserved: 



2 Substitute the true u(x, t) into the Lax-Friedrichs method (21) and use u t = cu x 
and u tt = c 2 u xx to find the coefficient of the numerical dissipation u xx . 

3 The difference equation Uj >n+ i = a mUj+m,n has growth factor G = a m e tmkAx 
Show consistency with e ickAt (first-order accuracy at least) when a m = 1 and 
^]mo m = cAt/Ax = r. 

4 The condition for second-order accuracy is ^m 2 a m = r 2 , from the Taylor se- 
ries. Check this for Lax- Wendroff with a$ = 1— r 2 ,ai = |(r 2 +r), a_i = \{j' 2 — r). 
With nonnegative coefficients a m > 0, the Schwarz inequality mv^V^') 2 — 
(J2m 2 a m )(J2a m ) becomes an equality r 2 = r 2 . This equality only happens if 
iTT'y/d^ = (constant) y/a^. Second-order is impossible with a m > 0, unless the 
difference equation has only one term Uj jn+ i = Uj +m , n . 

5 The Lax- Wendroff method has G = 1 — r 2 + r 2 cos kAx + ir sin kx. Square the 
real and imaginary parts to get (eventually!) \G\ 2 = 1 — (r 2 — r 4 )(l —cos kAx) 2 . 
Prove stability, that \G\ 2 < 1 if r 2 < 1. 

6 Suppose the coefficients in a linear differential equation change as t changes. 
The one-step solution operators become Sk and Rk, for the step from k At to 
(k + l)At. After n steps, products replace the powers S n and R n in U and u: 

U{n At) = S n -iS n ^2 ■ ■ ■ SiS w(0) and u{n At) = -R n _ii? n _ 2 . . . R\Rou(0) . 

Change the telescoping formula (25) to produce this U — u. Which parts are 
controlled by stability? Which parts by well-posedness (= stability of the dif- 
ferential equation) ? Consistency still controls Sk — Rk- 



T 



Problem Set 5.2 




and 




stay constant in time. 
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Even an unstable method will converge to the true solution u = e e for 
each separate frequency k. Consistency assures that the single-step growth 
factor G is 1 + ick At + O(At) 2 . Then for t = n At, 



How can we have convergence for each u(0) = e lhx and still prove divergence for 
a combination of frequencies u(0) = X^oo ° k e%kx ^ 

The upwind method with r > 1 is unstable because the CFL condition fails. 
By Problem 3, it does converge to e lk ( x+ct ) based on values of u(x, 0) = e %kx 
that do not reach as far as x + ct. The method must be finding a "correct" 
extrapolation of e tkx . So propose an initial u(x,0) for which convergence fails. 




which is convergence. 



